Seasonality in Lake Champlain Copepod Thermal Limits

### To Do 

Copepod Collection

Copepods were collected at approximately weekly intervals from Lake Champlain (Burlington Fishing Pier). Plankton was collected from the top 3 meters using a 250 um mesh net.

# Lake Champlain near Burlington, VT
siteNumber = "04294500"
ChamplainInfo = readNWISsite(siteNumber)
parameterCd = "00010"
startDate = "2023-01-01"
endDate = ""
#statCd = c("00001", "00002","00003", "00011") # 1 - max, 2 - min, 3 = mean

# Constructs the URL for the data wanted then downloads the data
url = constructNWISURL(siteNumbers = siteNumber, parameterCd = parameterCd, 
                       startDate = startDate, endDate = endDate, service = "uv")

temp_data = importWaterML1(url, asDateTime = T) %>% 
  mutate("date" = as.Date(dateTime)) %>% 
  select(date, "temp" = X_00010_00000)

Collections began in late May 2023. Several gaps are present, but collections have continued at roughly weekly intervals since then. Copepods from 39 collections were used to make a total of 1086 thermal limit measurements. Over this time period, collection temperatures ranged from 2.5 to 26.5°C.

There is substantial variation in thermal limits across the species collected. There is also some degree of variation within the species, with thermal limits increasing slightly during the summer.

## Daily values for the period examined by dataset
collection_conditions = temp_data %>%
  ungroup() %>% 
  group_by(date) %>% 
  summarise(mean_temp = mean(temp),
            med_temp = median(temp),
            var_temp = var(temp), 
            min_temp = min(temp), 
            max_temp = max(temp)) %>% 
  mutate("range_temp" = max_temp - min_temp,
         date = as.Date(date)) %>% 
  ungroup() %>%  
  filter(date >= (min(as.Date(full_data$collection_date)) - 7))

## Mean female thermal limits for each species, grouped by collection
species_summaries = full_data %>%  
  #filter(sex == "female") %>% 
  group_by(sp_name, collection_date, collection_temp) %>%  
  summarise("mean_ctmax" = mean(ctmax),
            "sample_size" = n(),
            "ctmax_st_err" = (sd(ctmax) / sqrt(sample_size)),
            "ctmax_var" = var(ctmax), 
            "mean_size" = mean(size),
            "size_st_err" = (sd(size) / sqrt(sample_size)),
            "size_var" = var(size)) %>%  
  ungroup() %>% 
  complete(sp_name, collection_date) %>% 
  arrange(desc(sample_size))

adult_summaries = full_data %>%  
  filter(sex == "female") %>% 
  group_by(sp_name, collection_date, collection_temp) %>%  
  summarise("mean_ctmax" = mean(ctmax),
            "sample_size" = n(),
            "ctmax_st_err" = (sd(ctmax) / sqrt(sample_size)),
            "ctmax_var" = var(ctmax), 
            "mean_size" = mean(size),
            "size_st_err" = (sd(size) / sqrt(sample_size)),
            "size_var" = var(size)) %>%  
  ungroup() %>% 
  complete(sp_name, collection_date) %>% 
  arrange(desc(sample_size))


ggplot() + 
  geom_vline(data = unique(select(full_data, collection_date)), 
             aes(xintercept = as.Date(collection_date)),
             colour = "grey90",
             linewidth = 1) + 
  geom_line(data = collection_conditions, 
            aes(x = as.Date(date), y = mean_temp),
            colour = "black", 
            linewidth = 2) + 
  # geom_errorbar(data = species_summaries,
  #               aes(x = as.Date(collection_date),
  #                   ymin = mean_ctmax - ctmax_st_err, ymax = mean_ctmax + ctmax_st_err,
  #                   colour = sp_name),
  #               position = position_dodge(width = 1),
  #               width = 5, linewidth = 1) +
  # geom_point(data = adult_summaries, 
  #            aes(x = as.Date(collection_date), y = mean_ctmax, colour = sp_name, size = sample_size)) + 
  geom_point(data = full_data, 
             aes(x = as.Date(collection_date), y = ctmax, colour = sp_name),
             size = 2, position = position_jitter(width = 1, height = 0)) + 
  scale_colour_manual(values = species_cols) + 
  labs(x = "Date", 
       y = "Temperature (°C)", 
       colour = "Species",
       size = "Sample Size") + 
  theme_matt() + 
  theme(legend.position = "right")

Size also varied, but primarily between rather than within species.

ggplot() + 
  geom_vline(data = unique(select(full_data, collection_date)), 
             aes(xintercept = as.Date(collection_date)),
             colour = "grey90",
             linewidth = 1) + 
  geom_line(data = collection_conditions, 
            aes(x = as.Date(date), y = mean_temp),
            colour = "black", 
            linewidth = 2) + 
  # geom_errorbar(data = species_summaries,
  #               aes(x = as.Date(collection_date), 
  #                   ymin = mean_ctmax - ctmax_st_err, ymax = mean_ctmax + ctmax_st_err,
  #                   colour = sp_name),
  #               position = position_dodge(width = 1),
  #               width = 5, linewidth = 1) + 
  geom_point(data = adult_summaries, 
             aes(x = as.Date(collection_date), y = mean_size * 40, colour = sp_name, size = sample_size),
             position = position_dodge(width = 1)) + 
  scale_colour_manual(values = species_cols) + 
  scale_y_continuous(
    name = "Temperature", # Features of the first axis
    sec.axis = sec_axis(~./40, name="Prosome Length (mm)"), # Add a second axis and specify its features
    breaks = c(0,5,10,15,20,25,30)
  ) + 
  labs(x = "Date", 
       y = "Temperature (°C)", 
       colour = "Species") + 
  theme_matt() + 
  theme(legend.position = "right")

Shown below is CTmax and body size for the species with the most data (Skistodiaptomus, L. minutus, L. sicilis, and Epischura), plotted against the day of the year for each sex/stage separately.

ctmax_feature = full_data %>%  
  mutate(doy = yday(collection_date)) %>% 
  filter(sp_name %in% c("Skistodiaptomus oregonensis", "Leptodiaptomus minutus", "Leptodiaptomus sicilis", "Epischura lacustris")) %>% 
  ggplot(aes(x = as.Date(collection_date), y = ctmax, colour = sp_name)) + 
  facet_grid(sp_name~sex) + 
  geom_point() + 
  scale_colour_manual(values = species_cols) + 
  labs(x = "Day of the Year", 
       y = "CTmax (°C)") + 
  theme_matt_facets() +
  theme(axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        legend.position = "none")

size_feature = full_data %>%  
  mutate(doy = yday(collection_date)) %>% 
  filter(sp_name %in% c("Skistodiaptomus oregonensis", "Leptodiaptomus minutus", "Leptodiaptomus sicilis", "Epischura lacustris")) %>% 
  ggplot(aes(x = as.Date(collection_date), y = size, colour = sp_name)) + 
  facet_grid(sp_name~sex) + 
  geom_point() + 
  scale_colour_manual(values = species_cols) + 
  labs(x = "Day of the Year", 
       y = "Size (mm)") + 
  theme_matt_facets() +
  theme(axis.text.x = element_text(angle = 300, hjust = 0, vjust = 0.5),
        legend.position = "none")

ggarrange(ctmax_feature, size_feature, common.legend = T, legend = "none")

adult_summaries %>% 
  ungroup() %>% 
  mutate(collection_num = as.numeric(factor(collection_date))) %>% 
  group_by(collection_date) %>%  
  arrange(collection_date) %>% 
  select(sp_name, collection_date, collection_num, sample_size) %>% 
  mutate(sample_size = replace_na(sample_size, 0)) %>% 
  mutate(total = sum(sample_size),
         percentage = sample_size / total,
         collection_date = lubridate::as_date(collection_date)) %>% 
  ggplot(aes(x = collection_date, y = percentage, fill = sp_name)) + 
  geom_area() + 
  scale_fill_manual(values = species_cols) + 
  scale_y_continuous(breaks = c(0,1)) + 
  labs(x = "Collection Date", 
       y = "Proportion", 
       fill = "Species") + 
  theme_minimal(base_size = 20) + 
  theme(panel.grid = element_blank(),
        axis.ticks = element_line())

Temperature Variability

Lake Champlain is highly seasonal, with both average temperatures and temperature variability changing throughout the year. These patterns in the experienced thermal environment may drive the observed variation in copepod thermal limits. However, the time period affecting copepod thermal limits is unknown. Depending the on the duration of time considered, there are large changes in the experienced environment, in particular regarding the temperature range and variance. Consider for example three time periods: the day of collection, one week prior to collection, and four weeks prior to collection. While the overall pattern is similar, we can see that, unsurprisingly, considering longer periods of time results in larger ranges and slightly changes the pattern of variance experienced.

## Daily values
daily_temp_data = temp_data %>%
  ungroup() %>% 
  group_by(date) %>% 
  summarise(mean_temp = mean(temp),
            med_temp = median(temp),
            var_temp = var(temp), 
            min_temp = min(temp), 
            max_temp = max(temp)) %>% 
  mutate("range_temp" = max_temp - min_temp)

day_prior_temp_data = temp_data %>% 
  ungroup() %>% 
  group_by(date) %>% 
  summarise(mean_temp = mean(temp),
            med_temp = median(temp),
            var_temp = var(temp), 
            min_temp = min(temp), 
            max_temp = max(temp)) %>% 
  mutate(date = date + 1) %>% 
  rename_with(.fn = ~ paste0("prior_day_", .x), .cols = c(-date))

daily_plot = daily_temp_data %>% 
  pivot_longer(cols = c(-date),
               names_to = "parameter", 
               values_to = "temp") %>% 
  ggplot(aes(x = date, y = temp, colour = parameter)) + 
  geom_line(linewidth = 1) + 
  scale_colour_manual(values = c(
    "mean_temp" = "olivedrab3",
    "med_temp" = "seagreen3",
    "max_temp" = "tomato",  
    "min_temp" = "dodgerblue",
    "range_temp" = "goldenrod3",
    "var_temp" = "darkgoldenrod1"
  )) + 
  scale_x_continuous(breaks = as.Date(c("2023-01-01", "2023-04-01", "2023-07-01"))) + 
  ggtitle("Daily Values") + 
  labs(y = "Temperature (°C)",
       x = "") + 
  theme_bw(base_size = 20) + 
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 270, hjust = 0, vjust = 0.5))

## Defining the function to get predictor values for periods of different lengths
get_predictors = function(daily_values, raw_temp, n_days){
  prefix = str_replace_all(xfun::numbers_to_words(n_days), pattern = " ", replacement = "-")
  
  mean_values = daily_values %>% 
    ungroup() %>% 
    mutate(mean_max = slide_vec(.x = max_temp, .f = mean, .before = n_days, .complete = T),
           mean_min = slide_vec(.x = min_temp, .f = mean, .before = n_days, .complete = T),
           mean_range = slide_vec(.x = range_temp, .f = mean, .before = n_days, .complete = T)) %>% 
    select(date, mean_max, mean_min, mean_range) %>% 
    rename_with( ~ paste(prefix, "day", .x, sep = "_"), .cols = c(-date))
  
  period_values = raw_temp %>% 
    mutate(mean = slide_index_mean(temp, i = date, before = days(n_days), 
                                   na_rm = T),
           max = slide_index_max(temp, i = date, before = days(n_days), 
                                 na_rm = T),
           min = slide_index_min(temp, i = date, before = days(n_days),
                                 na_rm = T),
           med = slide_index_dbl(temp, .i = date, .before = days(n_days), 
                                 na_rm = T, .f = median),
           var = slide_index_dbl(temp, .i = date, .before = days(n_days), 
                                 .f = var),
           range = max - min) %>%  
    select(-temp) %>%  
    distinct() %>% 
    rename_with( ~ paste(prefix, "day", .x, sep = "_"), .cols = c(-date))%>% 
    inner_join(mean_values, by = c("date")) %>%  
    drop_na()
  
  return(period_values)
}

# ## Getting predictor variables for different periods
# 
# ### Short (three days)
# three_day_temps = get_predictors(daily_values = daily_temp_data, 
#                                  raw_temp = temp_data, 
#                                  n_days = 3)
# 
# ### ONE WEEK
week_temps = get_predictors(daily_values = daily_temp_data,
                            raw_temp = temp_data,
                            n_days = 7)

week_plot = week_temps %>%
  pivot_longer(cols = c(-date),
               names_to = "parameter",
               values_to = "temp") %>%
  filter(parameter %in% c("seven_day_mean",
                          "seven_day_med",
                          "seven_day_max",
                          "seven_day_min",
                          "seven_day_var",
                          "seven_day_range")) %>%
  mutate(parameter = paste(word(parameter, start = 3, sep = fixed("_")), "_temp", sep = "")) %>%
  ggplot(aes(x = date, y = temp, colour = parameter)) +
  geom_line(linewidth = 1) +
  scale_colour_manual(values = c(
    "mean_temp" = "olivedrab3",
    "med_temp" = "seagreen3",
    "max_temp" = "tomato",
    "min_temp" = "dodgerblue",
    "range_temp" = "goldenrod3",
    "var_temp" = "darkgoldenrod1"
  )) +
  scale_x_continuous(breaks = as.Date(c("2023-01-01", "2023-04-01", "2023-07-01"))) +
  ggtitle("One Week") +
  labs(y = "Temperature (°C)",
       x = "") +
  theme_bw(base_size = 20) +
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 270, hjust = 0, vjust = 0.5))
# 
# 
# ### TWO WEEKS
# two_week_temps = get_predictors(daily_values = daily_temp_data, 
#                                 raw_temp = temp_data, 
#                                 n_days = 14)
# 
# two_week_plot = two_week_temps %>% 
#   pivot_longer(cols = c(-date),
#                names_to = "parameter", 
#                values_to = "temp") %>% 
#   filter(parameter %in% c("fourteen_day_mean",
#                           "fourteen_day_med",
#                           "fourteen_day_max", 
#                           "fourteen_day_min", 
#                           "fourteen_day_var",
#                           "fourteen_day_range")) %>% 
#   mutate(parameter = paste(word(parameter, start = 3, sep = fixed("_")), "_temp", sep = "")) %>% 
#   ggplot(aes(x = date, y = temp, colour = parameter)) + 
#   geom_line(linewidth = 1) + 
#   scale_colour_manual(values = c(
#     "mean_temp" = "olivedrab3",
#     "med_temp" = "seagreen3",
#     "max_temp" = "tomato",  
#     "min_temp" = "dodgerblue",
#     "range_temp" = "goldenrod3",
#     "var_temp" = "darkgoldenrod1"
#   )) + 
#   scale_x_continuous(breaks = as.Date(c("2023-01-01", "2023-04-01", "2023-07-01"))) + 
#   ggtitle("Two Weeks") + 
#   labs(y = "Temperature (°C)",
#        x = "") + 
#   theme_bw(base_size = 20) + 
#   theme(panel.grid = element_blank(),
#         axis.text.x = element_text(angle = 270, hjust = 0, vjust = 0.5))
# 
# 
# ### FOUR WEEKS
four_week_temps = get_predictors(daily_values = daily_temp_data,
                                 raw_temp = temp_data,
                                 n_days = 28)

four_week_plot = four_week_temps %>%
  pivot_longer(cols = c(-date),
               names_to = "parameter",
               values_to = "temp") %>%
  filter(parameter %in% c("twenty-eight_day_mean",
                          "twenty-eight_day_med",
                          "twenty-eight_day_max",
                          "twenty-eight_day_min",
                          "twenty-eight_day_var",
                          "twenty-eight_day_range")) %>%
  mutate(parameter = paste(word(parameter, start = 3, sep = fixed("_")), "_temp", sep = "")) %>%
  ggplot(aes(x = date, y = temp, colour = parameter)) +
  geom_line(linewidth = 1) +
  scale_colour_manual(values = c(
    "mean_temp" = "olivedrab3",
    "med_temp" = "seagreen3",
    "max_temp" = "tomato",
    "min_temp" = "dodgerblue",
    "range_temp" = "goldenrod3",
    "var_temp" = "darkgoldenrod1"
  )) +
  scale_x_continuous(breaks = as.Date(c("2023-01-01", "2023-04-01", "2023-07-01"))) +
  ggtitle("Four Weeks") +
  labs(y = "Temperature (°C)",
       x = "") +
  theme_bw(base_size = 20) +
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 270, hjust = 0, vjust = 0.5))
# 
# 
# ### EIGHT WEEKS
# eight_week_temps = get_predictors(daily_values = daily_temp_data, 
#                                   raw_temp = temp_data, 
#                                   n_days = 56)
# 
# eight_week_plot = eight_week_temps %>% 
#   pivot_longer(cols = c(-date),
#                names_to = "parameter", 
#                values_to = "temp") %>% 
#   filter(parameter %in% c("fifty-six_day_mean",
#                           "fifty-six_day_med",
#                           "fifty-six_day_max", 
#                           "fifty-six_day_min", 
#                           "fifty-six_day_var",
#                           "fifty-six_day_range")) %>% 
#   mutate(parameter = paste(word(parameter, start = 3, sep = fixed("_")), "_temp", sep = "")) %>% 
#   ggplot(aes(x = date, y = temp, colour = parameter)) + 
#   geom_line(linewidth = 1) + 
#   scale_colour_manual(values = c(
#     "mean_temp" = "olivedrab3",
#     "med_temp" = "seagreen3",
#     "max_temp" = "tomato",  
#     "min_temp" = "dodgerblue",
#     "range_temp" = "goldenrod3",
#     "var_temp" = "darkgoldenrod1"
#   )) + 
#   scale_x_continuous(breaks = as.Date(c("2023-01-01", "2023-04-01", "2023-07-01"))) + 
#   ggtitle("Eight Weeks") + 
#   labs(y = "Temperature (°C)",
#        x = "") + 
#   theme_bw(base_size = 20) + 
#   theme(panel.grid = element_blank(),
#         axis.text.x = element_text(angle = 270, hjust = 0, vjust = 0.5))
# 
ggarrange(daily_plot, week_plot, four_week_plot, 
          common.legend = T, nrow = 1, legend = "bottom")

The different time periods examined by this climate data highlights that the relationship between minimum and maximum temperatures changes based on the window examined. For example, minimum and maximum temperatures experienced over weekly intervals are closely linked, whereas there is a distinct seasonal cycle in the relationship between minimum and maximum temperatures experienced over periods of four weeks.

one_week_doy_data = week_temps %>% 
  mutate(doy = yday(date))

one_week_temp_circle = ggplot(one_week_doy_data, aes(x = seven_day_mean_max, y = seven_day_mean_min, colour = doy)) + 
  geom_point() + 
  scale_colour_gradient2(
    high = "dodgerblue4",
    mid = "coral2",
    low = "dodgerblue4",
    midpoint = 182.5) + 
  labs(x = "Max. Temp. (°C)",
       y = "Min. Temp. (°C)") + 
  labs(x = "Max. Temp. (°C)",
       y = "Min. Temp. (°C)") + 
  ggtitle("One Week") + 
  theme_matt()

four_week_doy_data = four_week_temps %>% 
  mutate(doy = yday(date))

four_week_temp_circle = ggplot(four_week_doy_data, aes(x = `twenty-eight_day_max`, y = `twenty-eight_day_min`, colour = doy)) + 
  geom_point() + 
  scale_colour_gradient2(
    high = "dodgerblue4",
    mid = "coral2",
    low = "dodgerblue4",
    midpoint = 182.5) + 
  labs(x = "Max. Temp. (°C)",
       y = "Min. Temp. (°C)") + 
  ggtitle("Four Week") + 
  theme_matt()

ggarrange(one_week_temp_circle, four_week_temp_circle,
          common.legend = T, legend = "bottom")

The thermal environment over any period of time may drive patterns in thermal acclimation. To explore the potential effects of different acclimation windows, we examined the correlation between thermal limits and different representations of the thermal environment for different periods of time. Shown below are the correlation coefficients for these relationships. Each facet shows the relationship for a different dimension of the thermal environment. Correlation coefficients are plotted for different durations, for species that were collected more than five times. Only data for mature female copepods was included.

We can see that, in general, copepods are responding to proximate cues from the thermal environment, with correlations generally dropping off substantially as acclimation window duration increases. An exception is Epischura lacustris, which appears to be responding to maximum temperatures experienced over a 20 day time period.

### Pulling predictors and measuring correlations for much finer timescales; 1-56 days

num_colls = full_data %>% 
  filter(sex == "female") %>% 
  select(collection_date, sp_name) %>%  
  distinct() %>%  
  count(sp_name) %>% 
  filter(n >= 5)

corr_vals = data.frame()

dur_vals = c(1:50)
for(i in dur_vals){
  
  duration_temps = get_predictors(daily_values = daily_temp_data, 
                                  raw_temp = temp_data, 
                                  n_days = i) %>% 
    filter(date %in% as_date(unique(full_data$collection_date)))
  
  corr_data = full_data %>%
    filter(sp_name %in% num_colls$sp_name) %>% 
    filter(sex == "female") %>% 
    mutate(collection_date = as.Date(collection_date)) %>% 
    inner_join(duration_temps, join_by(collection_date == date)) %>% 
    pivot_longer(cols = c(collection_temp, contains("day_")),
                 values_to = "value", 
                 names_to = "predictor") %>%  
    group_by(sp_name, predictor) %>% 
    summarise(correlation = cor.test(ctmax, value)$estimate,
              p.value = cor.test(ctmax, value)$p.value,
              ci_low = cor.test(ctmax, value)$conf.int[1],
              ci_high = cor.test(ctmax, value)$conf.int[2],
              .groups = "keep") %>% 
    filter(predictor != "collection_temp") %>% 
    mutate(sig = ifelse(p.value <0.05, "Sig.", "Non Sig.")) %>% 
    separate(predictor, "_day_", into = c(NA, "parameter")) %>% 
    mutate(duration = i)
  
  corr_vals = bind_rows(corr_vals, corr_data)
}

coll_corr = full_data %>%
  filter(sp_name %in% num_colls$sp_name) %>% 
  filter(sex == "female") %>% 
  group_by(sp_name) %>% 
  summarise(correlation = cor.test(ctmax, collection_temp)$estimate,
            p.value = cor.test(ctmax, collection_temp)$p.value,
            ci_low = cor.test(ctmax, collection_temp)$conf.int[1],
            ci_high = cor.test(ctmax, collection_temp)$conf.int[2]) %>% 
  mutate(sig = ifelse(p.value <0.05, "Sig.", "Non Sig.")) %>% 
  mutate(duration = 0,
         parameter = "coll_temp")

corr_vals = corr_vals %>%  
  mutate(duration = as.numeric(duration)) %>% 
  bind_rows(coll_corr)

corr_vals %>% 
  mutate(parameter = fct_relevel(parameter, c("min", "max", "range",
                                              "mean", "med", "var",
                                              "mean_min", "mean_max", "mean_range"))) %>% 
  ggplot(aes(x = duration, y = correlation, colour = sp_name)) + 
  facet_wrap(.~parameter) + 
  geom_hline(yintercept = 0) + 
  geom_point(size = 0.9) + 
  geom_line(linewidth = 1.5) + 
  scale_colour_manual(values = species_cols) + 
  labs(x = "Duration (days)",
       y = "Correlation", 
       colour = "Species") + 
  theme_matt_facets()

Shown here are the top three factors for each species.

# 
# corr_vals = full_data %>%
#   filter(sp_name %in% num_colls$sp_name) %>% 
#   filter(sex == "female") %>% 
#   mutate(collection_date = as.Date(collection_date)) %>% 
#   full_join(temp_predictors, join_by(collection_date == date)) %>% 
#   pivot_longer(cols = c(collection_temp, mean_temp:tail(names(.), 1)),
#                values_to = "value", 
#                names_to = "predictor") %>%  
#   group_by(sp_name, predictor) %>% 
#   summarise(correlation = cor.test(ctmax, value)$estimate,
#             p.value = cor.test(ctmax, value)$p.value,
#             ci_low = cor.test(ctmax, value)$conf.int[1],
#             ci_high = cor.test(ctmax, value)$conf.int[2]) %>% 
#   mutate(sig = ifelse(p.value <0.05, "Sig.", "Non Sig."))

corr_vals %>%  
  filter(sig == "Sig.") %>% 
  drop_na(correlation) %>% 
  group_by(sp_name) %>%
  arrange(desc(correlation)) %>% 
  slice_head(n = 3) %>% 
  select("Species" = sp_name, "Predictor" = parameter, "Duration" = duration, "Correlation" = correlation, "P-Value" = p.value) %>% 
  knitr::kable(align = "c")

Species	Predictor	Duration	Correlation	P-Value
Epischura lacustris	max	20	0.8926416	0.0000000
Epischura lacustris	max	19	0.8906963	0.0000000
Epischura lacustris	max	21	0.8874924	0.0000000
Leptodiaptomus minutus	max	6	0.7547703	0.0000000
Leptodiaptomus minutus	max	8	0.7547371	0.0000000
Leptodiaptomus minutus	max	7	0.7545529	0.0000000
Leptodiaptomus sicilis	range	23	0.3468229	0.0000000
Leptodiaptomus sicilis	range	24	0.3457819	0.0000000
Leptodiaptomus sicilis	var	21	0.3435719	0.0000000
Limnocalanus macrurus	max	8	0.5542854	0.0001397
Limnocalanus macrurus	max	7	0.5539836	0.0001412
Limnocalanus macrurus	mean_min	6	0.5530991	0.0001454
Senecella calanoides	mean_max	5	0.5434258	0.0445959
Senecella calanoides	max	4	0.5402865	0.0460906
Senecella calanoides	max	5	0.5402865	0.0460906
Skistodiaptomus oregonensis	max	2	0.7964802	0.0000000
Skistodiaptomus oregonensis	max	1	0.7903638	0.0000000
Skistodiaptomus oregonensis	mean_max	2	0.7896567	0.0000000

Phenotypic variation (like acclimation of thermal limits) is a physiological process. depending on the mechanistic underpinnings (changes in HSP expression, etc.), the amount of time it takes for an individual to acclimate may vary based on body size (larger species, more cells, more time required to acclimate). Shown here is the duration of the environmental acclimation window the copepods appear to be responding to.

mean_sizes = full_data %>% 
  filter(sex == "female") %>% 
  group_by(sp_name) %>%  
  summarise(mean_size = mean(size, na.rm = T))

corr_vals %>% 
  group_by(sp_name) %>% 
  filter(correlation == max(correlation)) %>%  
  inner_join(mean_sizes, by = "sp_name") %>% 
  select(sp_name, duration, mean_size) %>%  
  ggplot(aes(x = mean_size, y = duration)) + 
  geom_point(aes(colour = sp_name), 
             size = 4) + 
  scale_colour_manual(values = species_cols) + 
  labs(x = "Mean Female Size (mm)",
       y = "Acclimation Duration",
       colour = "Species") + 
  theme_matt() + 
  theme(legend.position = "right")

Trait Variation

# ctmax_plot = full_data %>%
#   mutate( #sp_name = str_replace(sp_name, pattern = " ",
#     #                              replacement = "\n"),
#     sp_name = fct_reorder(sp_name, ctmax, mean)) %>%
#   ggplot(aes(y = sp_name, x = ctmax)) +
#   geom_point(aes(colour= sp_name_sub),
#              position = position_dodge(width = 0.3),
#              size = 4) +
#   scale_colour_manual(values = species_cols) +
#   xlab(NULL) +
#   labs(y = "",
#        x = "CTmax (°C)",
#        colour = "Group") +
#   theme_matt() +
#   theme(legend.position = "none")
# 
# size_plot = full_data %>%
#   mutate(sp_name = fct_reorder(sp_name, ctmax, mean)) %>%
#   ggplot(aes(y = sp_name, x = size)) +
#   geom_point(aes(colour= sp_name_sub),
#              position = position_dodge(width = 0.3),
#              size = 4) +
#   scale_colour_manual(values = species_cols) +
#   labs(x = "Prosome Length (mm)",
#        y = "",
#        colour = "Group") +
#   guides(color = guide_legend(ncol = 1)) +
#   theme_matt(base_size = ) +
#   theme(legend.position = "right",
#         axis.text.y = element_blank(),
#         plot.margin = margin(0, 0, 0, 0,"cm"))
# 
# trait_plot = ctmax_plot + size_plot
# trait_plot

Shown below are the clutch size distributions for the three diaptomiid species, which produce egg sacs that allow for easy quantification of fecundity.

full_data %>%  
  drop_na(fecundity) %>%  
  ggplot(aes(x = fecundity, fill = sp_name_sub)) + 
  facet_wrap(.~sp_name_sub, ncol = 1) + 
  geom_histogram(binwidth = 2) + 
  scale_fill_manual(values = species_cols) + 
  labs(x = "Fecundity (# Eggs)") +
  theme_matt_facets() + 
  theme(legend.position = "none")

One of the main aims of this project is to examine the patterns and processes driving variation in upper thermal limits across these species of copepods.

Variation with temperature

We expect one of the primary drivers of copepod thermal limits to be temperature. The correlation analysis has shown that the copepods are generally (although not always) responding to the recent thermal environment. Shown below are thermal limits, body size, and fecundity values plotted against the temperature at the time of collection. Also shown is warming tolerance, calculated as the difference between upper thermal limit and the collection temperature.

We generally see an increase in thermal limits with increasing collection temperature, a slight decrease in body size, and variable relationships between temperature and fecundity. All species maintained some degree of buffer between environmental temperatures and upper thermal limits, but Epischura and L. minutus approached their upper thermal limits during the warmest collections during the summer.

ctmax_temp = ggplot(full_data, aes(x = collection_temp, y = ctmax, colour = sp_name)) + 
  geom_smooth(method = "lm", linewidth = 3) +
  geom_point(size = 3) + 
  labs(x = "Collection Temperature (°C)", 
       y = "CTmax (°C)",
       colour = "Species") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

size_temp = ggplot(filter(full_data, sex != "juvenile"), aes(x = collection_temp, y = size, colour = sp_name)) + 
  geom_smooth(method = "lm", linewidth = 3) +
  geom_point(size = 3) + 
  labs(x = "Collection Temperature (°C)", 
       y = "Length (mm)",
       colour = "Species")  + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

wt_temp = ggplot(full_data, aes(x = collection_temp, y = warming_tol, colour = sp_name)) + 
  geom_smooth(method = "lm", linewidth = 3) +
  geom_point(size = 3) + 
  labs(x = "Collection Temperature (°C)", 
       y = "Warming Tolerance (°C)",
       colour = "Species")  + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

eggs_temp = ggplot(full_data, aes(x = collection_temp, y = fecundity, colour = sp_name)) + 
  geom_smooth(method = "lm", linewidth = 3) +
  geom_point(size = 3) + 
  labs(x = "Collection Temperature (°C)", 
       y = "Fecundity (# Eggs)",
       colour = "Species")  + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

ggarrange(ctmax_temp, size_temp, wt_temp, eggs_temp, 
          common.legend = T, legend = "right")

full_data %>% 
  #filter(sex == "female") %>%  
  group_by(sp_name) %>% filter(n() > 5) %>% filter(!str_detect(sp_name, pattern = "kindti")) %>% 
  ggplot(aes(x = collection_temp, y = ctmax, colour = sp_name)) + 
  facet_wrap(sp_name~.) + 
  geom_point(size = 2, alpha = 0.8) + 
  geom_smooth(method = "lm", se = F, linewidth = 2) + 
  labs(x = "Collection Temp. (°C)", 
       y = "CTmax (°C)") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "none")

Temperature dependence is relatively weak in L. sicilis, especially at cooler temperatures. We will return to this feature later in the report, but for now we will note that there are two size morphs in this species, which appear to respond differently to decreases in temperature. There are significant differences between the morphs and how temperature affects CTmax.

morph_data = full_data %>% 
  filter(sex == "female" & sp_name == "Leptodiaptomus sicilis") %>% 
  mutate(morph = if_else(size > 0.89, "large", "small"))

ggplot(morph_data, aes(x = collection_temp, y = ctmax, colour = morph)) + 
  geom_point(size = 2, alpha = 0.8) + 
  geom_smooth(method = "lm", se = T, linewidth = 2) + 
  labs(x = "Collection Temp. (°C)", 
       y = "CTmax (°C)") + 
  theme_matt() + 
  theme(legend.position = "none")

morph.model = lm(data = morph_data, 
                 ctmax ~ collection_temp * morph)

knitr::kable(car::Anova(morph.model, type = "III", test = "F"))

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	11346.62713	1	4238.107970	0.0000000
collection_temp	115.44002	1	43.118298	0.0000000
morph	36.75982	1	13.730254	0.0002546
collection_temp:morph	14.76000	1	5.513046	0.0195762
Residuals	744.28551	278	NA	NA

#summary(morph.model)

#morph.em = emmeans::emmeans(morph.model, pairwise ~ morph)

#plot(morph.em)

Copepods spent several days in lab during experiments. Shown below are the CTmax residuals (taken from a model of CTmax against collection temperature) plotted against the time spent in lab before measurements were made. Individual regressions are shown for the residuals against days in lab for each collection. We can see clearly that thermal limits are fairly stable over time.

ggplot(ctmax_resids, aes(x = days_in_lab, y = resids, colour = sp_name, group = collection_date)) + 
  facet_wrap(sp_name~.) + 
  geom_point(size = 4, alpha = 0.5) + 
  geom_smooth(method = "lm", se = F, linewidth = 1) + 
  #scale_x_continuous(breaks = c(0:5)) + 
  labs(x = "Days in lab", 
       y = "CTmax Residuals") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt_facets() + 
  theme(legend.position = "none")

full.model = lme4::lmer(data = model_data,
                        ctmax ~ sex + temp_cent + size_cent +
                          (1 + days_in_lab + temp_cent + size_cent|sp_name:sex))

car::Anova(full.model)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ctmax
##             Chisq Df Pr(>Chisq)    
## sex        2.8465  2     0.2409    
## temp_cent 42.1916  1  8.275e-11 ***
## size_cent  2.0467  1     0.1525    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fixed = fixef(full.model)

model_coefs = coefficients(full.model)$`sp_name:sex` %>%  
  rownames_to_column(var = "species") %>% 
  separate(species, into = c("species", "sex"), sep = ":") %>% 
  select(species, sex, "intercept" = "(Intercept)", temp_cent, size_cent, days_in_lab)

ggplot(model_coefs, aes(x = intercept, y = temp_cent)) + 
  geom_smooth(method = "lm", colour = "black") +
  geom_point(aes(colour = species, shape = sex),
             size = 6) + 
  scale_colour_manual(values = species_cols) + 
  labs(x = "Species Intercept", 
       y = "ARR") +
  theme_matt() + 
  theme(legend.position = "right")

The term “acclimation response ratio” is often used to describe the effect of temperature on thermal limits. The ARR is calculated as the change in thermal limits per degree change in acclimation temperature. For our data, we will estimate ARR as the slope of CTmax against collection temperature. These slopes were taken from a regression of CTmax against collection temperature and body size. Two different model types were used, a simple linear regression and a mixed effects model. The estimated ARR values were generally highly similar between the model types used.

Sex and stage variation in thermal limits

Previous sections have generally lumped juvenile, female, and male individuals together. There may be important stage- or sex-specific differences in CTmax though. For several species, we have measurements for individuals in different stages or of different sexes.

sex_sample_sizes = ctmax_resids %>%  
  group_by(sp_name, sex) %>%  
  summarise(num = n()) %>%  
  pivot_wider(id_cols = sp_name,
              names_from = sex, 
              values_from = num,
              values_fill = 0) %>% 
  select("Species" = sp_name, "Juvenile" = juvenile, "Female" = female, "Male" = male)

knitr::kable(sex_sample_sizes, align = "c")

Species	Juvenile	Female	Male
Epischura lacustris	22	45	19
Leptodiaptomus minutus	10	209	33
Leptodiaptomus sicilis	31	282	78
Limnocalanus macrurus	2	42	36
Osphranticum labronectum	0	1	0
Senecella calanoides	1	14	8
Skistodiaptomus oregonensis	14	191	28

The female-male and female-juvenile comparisons show that there are generally no differences in thermal limits between these groups.

ctmax_resids %>% 
  filter(sp_name %in% filter(sex_sample_sizes, Male > 0, Female > 0)$Species & 
           sex != "juvenile") %>% 
  ggplot(aes(x = sex, y = resids, colour = sp_name, group = sp_name)) + 
  facet_wrap(sp_name~., ncol = 2) + 
  geom_smooth(method = "lm", se = F, linewidth = 1) + 
  geom_point(size = 3,
             alpha = 0.5,
             position = position_jitter(height = 0, width = 0.05)) +  
  labs(x = "Sex", 
       y = "CTmax Residuals") + 
  scale_colour_manual(values = species_cols) + 
  theme_bw(base_size = 18) + 
  theme(legend.position = "none", 
        panel.grid = element_blank())

ctmax_resids %>% 
  filter(sp_name %in% filter(sex_sample_sizes, Juvenile > 0 & Female > 0)$Species & 
           sex != "male") %>% 
  ggplot(aes(x = sex, y = resids, colour = sp_name, group = sp_name)) + 
  facet_wrap(sp_name~., ncol = 2) + 
  geom_smooth(method = "lm", se = F, linewidth = 1) + 
  geom_point(size = 3,
             alpha = 0.5,
             position = position_jitter(height = 0, width = 0.05)) +  
  labs(x = "Sex", 
       y = "CTmax (°C)") + 
  scale_colour_manual(values = species_cols) + 
  theme_bw(base_size = 18) + 
  theme(legend.position = "none", 
        panel.grid = element_blank())

Trait Correlations and Trade-offs

A relationship between size and upper thermal limits has been suggested in a wide range of other taxa. Shown below are the measured upper thermal limits plotted against prosome length. The overall relationship (inclusive of all species) is shown as the black line in the background. Regressions for each individual species are also shown. Across the entire assemblage, there is a strong decrease in thermal limits with increasing size.

full_data %>% 
  #filter(sex == "female") %>%  
  ggplot( aes(x = size, y = ctmax, colour = sp_name)) + 
  geom_smooth(data = full_data, 
              aes(x = size, y = ctmax),
              method = "lm", 
              colour ="black", 
              linewidth = 2.5) + 
  geom_point(size = 2, alpha = 0.3) + 
  geom_smooth(method = "lm", se = F, linewidth = 2) + 
  labs(x = "Length (mm)", 
       y = "CTmax (°C)",
       colour = "Species") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

Shown here is the relationship for each species individually.

full_data %>% 
  #filter(sex == "female") %>%  
  group_by(sp_name) %>% filter(n() >2) %>% filter(!str_detect(sp_name, pattern = "kindti")) %>% 
  ggplot( aes(x = size, y = ctmax, colour = sp_name)) + 
  facet_wrap(sp_name~., scales = "free", nrow = 2) + 
  geom_point(size = 2, alpha = 0.8) + 
  geom_smooth(method = "lm", se = F, linewidth = 2) + 
  labs(x = "Length (mm)", 
       y = "CTmax (°C)",
       colour = "Species") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "none")

Shown below is the relationship between mean size and mean thermal limits for females of each species. We see that larger species within the community tend to have a lower thermal limit than smaller species.

full_data %>% 
  group_by(sp_name, sex) %>% 
  summarize(mean_ctmax = mean(ctmax, na.rm = T),
            mean_size = mean(size, na.rm = T)) %>% 
  #filter(sex == "female") %>% 
  ggplot(aes(x = mean_size, y = mean_ctmax)) + 
  geom_smooth(method = "lm", se = F, linewidth = 2, colour = "black") + 
  geom_point(aes(colour = sp_name, shape = sex),
             size = 5) + 
  labs(x = "Length (mm)", 
       y = "CTmax (°C)",
       colour = "Species") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

Shown here is the relationship between fecundity and size, showing the classic pattern of increasing egg production with increasing size.

ctmax_resids %>%  
  drop_na(fecundity) %>% 
  ggplot(aes(x = size, y = fecundity, colour = sp_name)) + 
  geom_smooth(method = "lm", se = F, linewidth = 2) + 
  geom_point(size = 2, alpha = 0.5) + 
  labs(x = "Prosome length (mm)", 
       y = "Fecundity (# Eggs)",
       colour = "Species") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

Individuals may also allocate energy to different fitness related traits, prioritizing reproductive output over environmental tolerance, for example. Shown below is the relationship between CTmax residuals (again, controlling for the effects of collection temperature) against fecundity. We can see clearly that individuals with increased fecundity are not decreasing thermal limits, suggesting that there is no energetic trade-off between these traits.

ctmax_resids %>%  
  drop_na(fecundity) %>% 
  ggplot(aes(x = resids, y = fecundity, colour = sp_name)) + 
  geom_smooth(method = "lm", se = F, linewidth = 2) + 
  geom_point(size = 2, alpha = 0.5) + 
  labs(x = "CTmax Residuals", 
       y = "Fecundity (# Eggs)") + 
  scale_colour_manual(values = species_cols) + 
  theme_matt() + 
  theme(legend.position = "right")

fitness.model = lm(data = ctmax_resids, 
   fecundity ~ resids * sp_name)

car::Anova(fitness.model)
## Anova Table (Type II tests)
## 
## Response: fecundity
##                Sum Sq  Df  F value    Pr(>F)    
## resids            4.5   1   0.2756  0.600038    
## sp_name        8291.3   2 252.5993 < 2.2e-16 ***
## resids:sp_name  197.7   2   6.0234  0.002739 ** 
## Residuals      4693.8 286                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans::emtrends(fitness.model,  var = "resids","sp_name")
##  sp_name                     resids.trend    SE  df lower.CL upper.CL
##  Leptodiaptomus minutus             0.508 0.284 286  -0.0515    1.068
##  Leptodiaptomus sicilis            -0.129 0.250 286  -0.6217    0.364
##  Skistodiaptomus oregonensis       -1.213 0.407 286  -2.0147   -0.411
## 
## Confidence level used: 0.95

Other patterns in variation

Leptodiaptomus sicilis is the most abundant species during the winter. There was a large shift in the size of mature females towards the end of December. These large and small individuals are the same species (confirmed via COI sequencing), suggesting this shift may reflect a transition from one generation to another and that, unlike in many other lakes, there are two generations of L. sicilis per year in Lake Champlain. This size difference may be caused by differences in the developmental environments. For example, individuals developing in January grow up at very low temperatures, and therefore may reach larger sizes. These individuals oversummer in deep waters, then re-emerge in October and produce a new generation. Water temperatures are still fairly high through November, which results in a generation of smaller individuals, which mature in time to produce a new generation in January.

Shown below is the distribution of pairwise distances between COI sequences of large and small morphs. Distances in both within- and across-morph comparisons are small.

ind_dist = ape::dist.dna(sic_dnabin, model = "raw") %>% as.matrix %>% 
  as_tibble() %>%
  mutate("ind1" = colnames(.)) %>% 
  pivot_longer(-ind1, names_to = "ind2", values_to = "dist") %>%
  mutate(ind1 = factor(ind1),
         ind2 = factor(ind2)) %>% 
  filter(!(ind1 == "sore1" | ind2 == "sore1")) %>% 
  mutate(
    ind1 = case_when(
      ind1 == "S1" ~ "small1",
      ind1 == "S3" ~ "small3",
      ind1 == "lsic3" ~ "small4",
      ind1 == "lsic5" ~ "small6",
      ind1 == "lsic9" ~ "small8",
      ind1 == "lsic10" ~ "small9",
      ind1 == "lsic11" ~ "small10",
      ind1 == "L1" ~ "large1",
      ind1 == "L2" ~ "large2",
      ind1 == "L3" ~ "large3",
      ind1 == "lsic1" ~ "large4",
      ind1 == "lsic2" ~ "large5",
      ind1 == "lsic7" ~ "large6",
      ind1 == "lsic8" ~ "large7"),
    ind2 = case_when(
      ind2 == "S1" ~ "small1",
      ind2 == "S3" ~ "small3",
      ind2 == "lsic3" ~ "small4",
      ind2 == "lsic5" ~ "small6",
      ind2 == "lsic9" ~ "small8",
      ind2 == "lsic10" ~ "small9",
      ind2 == "lsic11" ~ "small10",
      ind2 == "L1" ~ "large1",
      ind2 == "L2" ~ "large2",
      ind2 == "L3" ~ "large3",
      ind2 == "lsic1" ~ "large4",
      ind2 == "lsic2" ~ "large5",
      ind2 == "lsic7" ~ "large6",
      ind2 == "lsic8" ~ "large7"),
    'comparison' = case_when(
      str_detect(ind1, pattern = "large") & str_detect(ind2, pattern = "large") ~ "within",
      str_detect(ind1, pattern = "small") & str_detect(ind2, pattern = "small") ~ "within", 
      str_detect(ind1, pattern = "large") & str_detect(ind2, pattern = "small") ~ "across",
      str_detect(ind1, pattern = "small") & str_detect(ind2, pattern = "large") ~ "across"
    )) 

ggplot(ind_dist, aes(dist, fill = comparison)) +
  geom_histogram(binwidth = 0.005) + 
  labs(x = "Distance") + 
  theme_matt()

full_data %>%  
  filter(sp_name == "Leptodiaptomus sicilis") %>% 
  filter(sex == "female") %>% 
  group_by(collection_date) %>% 
  mutate(size_center = scale(size, center = T, scale = F)) %>% 
  ggplot(aes(y = collection_date, x = size, fill = collection_temp)) + 
  geom_density_ridges(bandwidth = 0.04) + 
  geom_vline(xintercept = 0.89) + 
  labs(x = "Size (mm)",
       y = "Date", 
       fill = "Coll. Temp. (°C)") + 
  theme_matt() + 
  theme(legend.position = "right",
        axis.text.y = element_text(size = 12))

full_data %>%  
  filter(sp_name == "Leptodiaptomus minutus") %>% 
  filter(sex == "female") %>% 
  group_by(collection_date) %>% 
  mutate(size_center = scale(size, center = T, scale = F)) %>% 
  ggplot(aes(y = collection_date, x = size, fill = collection_temp)) + 
  geom_density_ridges(bandwidth = 0.04) + 
  geom_vline(xintercept = 0.69) + 
  labs(x = "Size (mm)",
       y = "Date", 
       fill = "Coll. Temp. (°C)") + 
  coord_cartesian(xlim = c(0.5,0.9)) + 
  theme_matt() + 
  theme(legend.position = "right",
        axis.text.y = element_text(size = 12))

if(predict_vuln == F){
  knitr::knit_exit()
}